Solution for Homework #2¶

IEDA/ELEC3180¶

TAs: Amir Javaheri and Yifan Yu¶

First ckeck if all necessary packages and modules (e.g., yfinance) are installed

In [2]:
# !pip install yfinance
In [3]:
import numpy as np
import cvxpy as cp

import pandas as pd
from pandas_datareader import data as pdr
import yfinance as yf
yf.pdr_override()

from matplotlib import animation
from IPython.display import HTML
from matplotlib import rc

from sklearn.preprocessing import scale

import seaborn as sns
import matplotlib.pyplot as plt
plt.style.use('ggplot')

Problem 1¶

a¶

In this problem we aim to model the log-prices of a set of stocks $\mathbf{x}_t$ as linear combination of factors $f_t$. As mentioned in the homework, we have a multi-factor model in which the S&P500 index and the Bitcoin price are used as factors. The coefficients $\boldsymbol{\beta}_i$ measures the degree of variability of the $i$-th stock with respect the factor $f_t$. The larger $\boldsymbol{\beta}_i$, the higher the sensitivity of the $i$-th stock will be, meaning that small changes in $f_t$ will lead to large changes in $\mathbf{x}_{t,i}$, which may be undesirable in many scenarios. $\mathbf{\alpha}_i$ is the part of the log-price of the $i$-th stock $\mathbf{x}_{t,i}$ which is not explained by the factor $f_t$. Hence, the higher $\mathbf{\alpha}_i$ the higher is the return of stock $i$ with respect to the factors $f_t$. Fianlly, $\boldsymbol{\epsilon}_i$ denotes the risk associated with the $i$-th stock, i.e., risks that cannot be explained by the factors and are usually known as unsystematic risk.

b¶

First, we import the data using the pdr.get_data_yahoomodule. We then merge the data corresponding to the stocks and the factors (using join method), interpolate the data and finally detach those two parts. Fianlly we compute the log returns via .apply(np.log).apply(np.diff) method.

In [4]:
stocks = ['TSLA', 'AMZN', 'EBAY', 'AAPL', 'MSFT', 'META', 'GOOGL',
          'NFLX', 'IBM', 'NVDA']
          
indices = ['^GSPC', 'BTC-USD']
prices = pdr.get_data_yahoo(stocks, start="2019-01-01", end="2022-12-31")[['Adj Close']]
factors = pdr.get_data_yahoo(indices, start="2019-01-01", end="2022-12-31")[['Adj Close']]

data_merged = factors.join(prices)
data_merged = data_merged.interpolate()

prices = data_merged['Adj Close'][stocks].fillna(1)
factors = data_merged['Adj Close'][indices].fillna(1)


log_ret_prices = prices.apply(np.log).apply(np.diff)
log_ret_factors = factors.apply(np.log).apply(np.diff)

[*********************100%***********************]  10 of 10 completed
[*********************100%***********************]  2 of 2 completed
In [5]:
print(prices.head())

                 TSLA       AMZN       EBAY       AAPL       MSFT        META  \
Date                                                                            
2019-01-01   1.000000   1.000000   1.000000   1.000000   1.000000    1.000000   
2019-01-02  20.674667  76.956497  27.043747  38.047050  96.632675  135.679993   
2019-01-03  20.024000  75.014000  26.546926  34.257282  93.077736  131.740005   
2019-01-04  21.179333  78.769501  27.156227  35.719700  97.406708  137.949997   
2019-01-05  21.563111  79.671501  27.312461  35.693200  97.448120  137.983332   

                GOOGL        NFLX        IBM       NVDA  
Date                                                     
2019-01-01   1.000000    1.000000   1.000000   1.000000  
2019-01-02  52.734001  267.660004  89.489304  33.799732  
2019-01-03  51.273499  271.200012  87.702805  31.757654  
2019-01-04  53.903500  297.570007  91.128258  33.792286  
2019-01-05  53.867667  303.493337  91.343155  34.388615

c¶

$\mathbf{X}\in \mathbb{R}^{N\times T}$ represents the log-returns of the stocks

$\mathbf{F} = [\mathbf{1}, \mathbf{f}^\top]^\top \in \mathbb{R}^{3\times T}$ where $\mathbf{f}\in \mathbb{R}^{2\times T}$ represents the factors log returns.

In [6]:
X = log_ret_prices.to_numpy().T
(N, T) = X.shape

F = np.vstack((np.ones((1,T)), log_ret_factors.to_numpy().T))
print((N,T))

(10, 1459)
In [7]:
Gamma = cp.Variable((N,3))
obj = cp.sum_squares(X - Gamma@ F)
problem = cp.Problem(cp.Minimize(obj))
result = problem.solve(verbose=False)

# alpha = Gamma.value[:,0]
# beta = Gamma.value[:,1:3]
# print('alpha = \n', alpha)
# print('beta = \n', beta)

Gamma = pd.DataFrame(Gamma.value, columns = ["alpha", "beta1","beta2"], index=stocks)
print(Gamma)

          alpha     beta1     beta2
TSLA   0.000937  0.388590  0.163798
AMZN  -0.000170  0.555907  0.061970
EBAY   0.000125  0.422265  0.035914
AAPL   0.000629  0.466464  0.066674
MSFT   0.000382  0.585149  0.061792
META  -0.000340  0.628715  0.067122
GOOGL  0.000138  0.507886  0.059752
NFLX  -0.000206  0.714680  0.060979
IBM    0.000124  0.574797  0.007824
NVDA   0.000717  0.452343  0.140429

d¶

Now using the closed form solution \begin{equation} \boldsymbol{\Gamma}^* = \underset{\boldsymbol{\Gamma}}{\text{argmin}}\, \|\mathbf{X} - \boldsymbol{\Gamma}\mathbf{F}\|_F^2 = \mathbf{X} \mathbf{F}^\top (\mathbf{F}\mathbf{F}^\top)^{-1} \end{equation}

In [8]:
Gamma = np.linalg.solve( F@F.T, F@X.T).T
# alpha = Gamma[:,0]
# beta = Gamma[:,1:3]

# print('closed-form alpha =\n ', alpha)
# print('closed-form beta =\n', beta)

Gamma = pd.DataFrame(Gamma, columns = ["alpha", "beta1","beta2"], index=stocks)
print('closed-form solution:')
print(Gamma)

closed-form solution:
          alpha     beta1     beta2
TSLA   0.000937  0.388590  0.163798
AMZN  -0.000170  0.555907  0.061970
EBAY   0.000125  0.422265  0.035914
AAPL   0.000629  0.466464  0.066674
MSFT   0.000382  0.585149  0.061792
META  -0.000340  0.628715  0.067122
GOOGL  0.000138  0.507886  0.059752
NFLX  -0.000206  0.714680  0.060979
IBM    0.000124  0.574797  0.007824
NVDA   0.000717  0.452343  0.140429
In [9]:
error = np.linalg.norm(X - Gamma@ F, 'fro')
print("Error closed form: ", error)

Error closed form:  2.0854391368107734

e¶

For this part we use fig.add_subplot(projection='3d') method to draw a 3D scatter plot

In [10]:
plt.rcParams.update(plt.rcParamsDefault)
In [11]:
fig = plt.figure()
ax = fig.add_subplot(projection='3d')

Gamma = Gamma.to_numpy()
for i in range(N):
    ax.scatter(Gamma[i,1],Gamma[i,2],Gamma[i,0],color='b')
    ax.text(Gamma[i,1],Gamma[i,2],Gamma[i,0],  '%s' % ( stocks[i]), size=10, zorder=1, color='k')


ax.set_xlabel('beta_1')
ax.set_ylabel('beta_2')
ax.set_zlabel('alpha')
plt.show()

We then use animation module to animate the figure

In [12]:
rc('animation', html='jshtml')
In [13]:
def animate(frame):
  ax.view_init(30, frame/4)
  plt.pause(.001)
  return fig
anim = animation.FuncAnimation(fig, animate, frames=200, interval=2)
# anim.save('animationBrownianMotion2d.gif', writer='pillow', fps=60)
anim
Out[13]:

f¶

The stocks with higher $\alpha$ and lower $\mathbf{\beta}$s would be less riskier to invest on. Hene, according to the above figure, stocks like 'TSLA', 'NVDA' and 'AAPL' have high $\alpha$s while also having small $\beta_1$s and $\beta_2$s, and it would probably be less riskier to invest on them.

Alternatively, one can compute the ratio of $\alpha_i/\|\mathbf{\beta}_i \|$ for each stock and choose the ones with the highest ratios as follows

In [14]:
alpha_over_beta = np.zeros((N,1))
for i in range(N):
  alpha_over_beta[i] = Gamma[i,0]/ np.linalg.norm(Gamma[i,1:3])

alpha_over_beta = pd.DataFrame(alpha_over_beta, columns = ["alpha/ beta"], index=stocks)
print( alpha_over_beta.sort_values(by=['alpha/ beta']) )

       alpha/ beta
META     -0.000538
AMZN     -0.000304
NFLX     -0.000288
IBM       0.000215
GOOGL     0.000269
EBAY      0.000294
MSFT      0.000650
AAPL      0.001336
NVDA      0.001513
TSLA      0.002222

e¶

We just use the code given in the Py-session_iid_modeling.html and we choose the number of factor $K=2$

In [15]:
K = 2
alpha = np.mean(X, axis=1)
X_sc = scale(X, with_std=False)
Sigma_prev = np.zeros((N, N))
Psi = np.zeros((N, N))
Sigma = np.cov(X_sc) 
while (np.linalg.norm(Sigma - Sigma_prev,'fro')/np.linalg.norm(Sigma,'fro') > 1e-5):
  values, vectors = np.linalg.eig(Sigma - Psi)
  B = vectors[:, :K] * np.sqrt(values[:K])
  Psi = np.diag(np.diag(Sigma - B @ B.T))
  Sigma_prev = Sigma
  Sigma = B @ B.T + Psi
In [16]:
Gamma = np.column_stack([alpha, B]).real
pd.DataFrame(Gamma, index=stocks, columns=[['alpha','beta_1','beta_2']])
Out[16]:
alpha beta_1 beta_2
TSLA 0.003299 -0.032772 0.017241
AMZN 0.003037 0.005894 -0.001157
EBAY 0.002549 -0.021648 -0.007147
AAPL 0.003335 -0.012684 -0.004288
MSFT 0.003754 0.010975 -0.002488
META 0.003283 0.021459 -0.000011
GOOGL 0.003071 -0.003735 -0.004319
NFLX 0.003898 0.039980 0.010717
IBM 0.003383 0.009151 -0.007620
NVDA 0.003416 -0.016620 -0.000927

Now, we can obtain the factors as follows \begin{equation} \mathbf{F} = \mathbf{\hat{\Lambda}}_1^{-1/2} \mathbf{\hat{\Gamma}}_1^\top (\mathbf{X} - \mathbf{\alpha} \mathbf{1}^\top) \end{equation}

In [17]:
B_inv = vectors[:, :K].real * np.sqrt(1/values[:K].real)
factor =  B_inv.T @ (X - alpha.reshape((10,1))@ np.ones((1,T)) )
factor.shape
Out[17]:
(2, 1459)
In [18]:
pd.DataFrame(factor.T, index=log_ret_prices.index, columns=[['factor1','factor2']])
Out[18]:
factor1 factor2
0 36.251165 3.273816
1 0.704806 1.088491
2 0.399514 1.541227
3 -0.034374 0.752570
4 -0.034361 0.736975
... ... ...
1454 0.204206 -1.023567
1455 0.212150 -1.053837
1456 -0.488302 1.200260
1457 -0.309674 1.940839
1458 0.044385 0.708391

1459 rows × 2 columns

In [19]:
F = np.vstack((np.ones((1,T)), factor))
error = np.linalg.norm(X - Gamma@ F, 'fro')
print("Error: ", error)

Error:  13.251833836956733

The previous model had less approximation error. So it seemed to be a better way of describing the stock prices via the linear model.

Problem 2¶

a¶

We first divide the data into training ans test sets.

In [20]:
prices_train = prices.loc['2019-01-01':'2021-12-31', :][['AAPL','MSFT','GOOGL']]
prices_test = prices.loc['2022-01-01':'2022-12-31', :][['AAPL','MSFT','GOOGL']]
In [21]:
log_returns_train = prices_train.apply(np.log).apply(np.diff)
log_returns_test = prices_test.apply(np.log).apply(np.diff)

Then, we define a function to design a portfolio based on the shrinkage estimator. The function takes the mean and covariance of the training set as well as the parameters $\rho_1$ and $\rho_2$ as input and gives the Sharpe ratio of the designed portfolio on the training data as the output result.

In [22]:
def design_portfolio(beta, mu, Sigma):
  w = cp.Variable(len(mu))
  variance = cp.quad_form(w, Sigma)
  exp_return = w @ mu
  problem = cp.Problem(cp.Minimize(variance), [exp_return >= beta, w >= 0, cp.sum(w) == 1])   
  result = problem.solve()
  return w.value

b¶

In this part, we design a porfolio based on the shrinkage estimator for each pair of $(\rho_1,\rho_2)$ (using the training data). In the following we first estimate the mean and sample covariance of the training data and then we use different pairs of $(\rho_1, \rho_2)$ to refine our estimations using the shrinkage method. Then the refined estimations of mean and covariance are used to design a portfolio.

In [23]:
rho = np.array([0, 0.2, 0.4, 0.6, 0.8, 1])
mu_train = np.mean(log_returns_train)
SCM_train = np.cov(log_returns_train.T)

/usr/local/lib/python3.9/dist-packages/numpy/core/fromnumeric.py:3472: FutureWarning: In a future version, DataFrame.mean(axis=None) will return a scalar mean over the entire DataFrame. To retain the old behavior, use 'frame.mean(axis=0)' or just 'frame.mean()'
  return mean(axis=axis, dtype=dtype, out=out, **kwargs)
In [24]:
def design_portfolio_with_shrinkage(beta, rho_1, rho_2, mu_hat, Sigma_hat):
    N = Sigma_hat.shape[0]
    Sigma = rho_1 * Sigma_hat + (1-rho_1) * np.diag(np.diag(Sigma_hat) )
    mu = rho_2 * mu_hat + (1-rho_2) * np.sum(mu_hat) * np.ones(N) / N
    w  = design_portfolio(beta, mu, Sigma)
    expected_return = np.dot(w, mu)
    volatility = np.sqrt(w.T @ Sigma @ w)
    sharpe_ratio = expected_return/volatility
    return sharpe_ratio, w

We choose $\beta=0.002$ in portfolio design

In [25]:
sharpe_ratios = []
portfolios = []
beta = 0.002
for rho_1 in rho:
    for rho_2 in rho:
        sr, w = design_portfolio_with_shrinkage(beta, rho_1, rho_2, mu_train, SCM_train)
        sharpe_ratios.append(sr)
        portfolios.append(w)
sharpe_ratios = np.array(sharpe_ratios).reshape(6,6)

Next, we find the pair which gives the highest Sharpe ratio

In [26]:
idx_max_sharpe_ratio = np.argmax(sharpe_ratios)
i, j = int(idx_max_sharpe_ratio / len(rho)), idx_max_sharpe_ratio % len(rho)
rho[i], rho[j]
Out[26]:
(0.0, 0.0)

Finally we plot the heatmap of the Sharpe ratios computed for different portfolios designed on the training data, via different pairs of $(\rho_1,\rho_2)$ in the shrinkage estimator.

In [27]:
sns.heatmap(sharpe_ratios, xticklabels=rho, yticklabels=rho)
plt.xlabel(r'$\rho_2$')
plt.ylabel(r'$\rho_1$')
plt.show()

c¶

In this part, we use the weights of the portfolio designed on the training data (previous part) for the pair $(\rho_1=0,\rho_2=0)$ (the pair which gave the highest Sharpe ratio) to compute the Sharpe ratio on the test data.

Note: We never design a portfolio based on the test data (since as the name suggests, this data is only used for test and is not available for training). We just use the mean and covariance of the test data to compute the Sharpe ratio as follows.

In [28]:
mu_test = np.mean(log_returns_test)
SCM_test = np.cov(log_returns_test.T)

The Sharpe ratio obtained on the test data, as well as the weights of the designed portfolio are as follows

In [29]:
w = portfolios[idx_max_sharpe_ratio]
expected_return = np.dot(w, mu_test)
volatility = np.sqrt(w.T @ SCM_test @ w)
test_sharpe_ratio = expected_return/volatility
print("SR: ", test_sharpe_ratio, "  w: ", w)

SR:  -0.06143629912066783   w:  [0.40189126 0.25729447 0.34081427]

d¶

Here, we again use the weights of the portfolio designed on the training data for the pair $(\rho_1=1,\rho_2=1)$ (no shrinkage) to compute the Sharpe ratio on the test set.

In [30]:
w = portfolios[-1]
expected_return = np.dot(w, mu_test)
volatility = np.sqrt(w.T @ SCM_test @ w)
test_sharpe_ratio_no_shrinkage = expected_return/volatility
print("SR: ", test_sharpe_ratio_no_shrinkage, "  w: ", w)

SR:  -0.04901758717964652   w:  [1.00000000e+00 4.33665699e-24 2.81887149e-24]

According to the results, using the shrinkage estimator in this case does not help since the Sharpe ratio is higher for the case of having no shrinkage estimator used.